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ABSTRACT 



Aims. This work presents a new hydrodynamical algorithm to study astrophysical detonations. A prime motivation of this development 
is the description of a carbon detonation in conditions relevant to superbursts, which are thought to result from the propagation of a 
detonation front around the surface of a neutron star in the carbon layer underlying the atmosphere. 

Methods. The algorithm we have developed is a finite- volume method inspired by the original MUSCL scheme of van Leer (1979). 
The algorithm is of second-order in the smooth part of the flow and avoids dimensional splitting. It is applied to some test cases, and 
the time-dependent results are compared to the corresponding steady state solution. 

Results. Our algorithm proves to be robust to test cases, and is considered to be reliably applicable to astrophysical detonations. The 
preliminary one-dimensional calculations we have performed demonstrate that the carbon detonation at the surface of a neutron star 
is a multiscale phenomenon. The length scale of liberation of energy is 10 6 times smaller than the total reaction length. We show that 
a multi-resolution approach can be used to solve all the reaction lengths. This result will be very useful in future multi-dimensional 
simulations. We present also thermodynamical and composition profiles after the passage of a detonation in a pure carbon or mixed 
carbon-iron layer, in thermodynamical conditions relevant to superbursts in pure helium accretor systems. 
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1. Introduction 



Methods: numerical - Shock waves - Stars: neutron - X-rays: bursts - Nuclear reactions, nucleosyn- 



Superbursts have been discovered by long term monitoring of 
the X-ray sky with instruments such as RXTE and BeppoSAX. 
Compared to normal type I X-ray bursts, they are 1000 times 
more energetic (integrated burst energies of about 10 42 ergs), 
1000 times longer (they last from hours to half a day), and have 
recurrence times of the order of years. They are very rare, only 
13 such events having been found from 8 sources (for reviews 
see Kuulkers 2004, Cumming 2005, and references therein). 

They also exhibit similarities with normal type I X-ray 
bursts, like a rapid rise in the light curve, a quasi-exponential 
decay, and a hardening of the spectrum during the rise followed 
by a softening during the decay, which are well represented by a 
blackbody model with an effective temperature growing during 
the rise and decreasing during the decay phase (Kuulkers 2004). 
This leads to the suggestion that superbursts, like normal type I 
bursts, are thermonuclear in origin (Cornelisse et al. 2000). The 
current view is that superbursts are due to the thermally unstable 
ignition of 12 C at densities of about 10 8 - 10 y g ctrT 3 (Cumming 
I2TJ0T1 Strohmayer & Brown 120021 Cumming|2005]l. 

Most superbursts have been observed in systems accreting a 
mix of H and He (Kuulkers |2004| l . In these H/He accretors a very 
small amount of 12 C remains after the combustion of H and He 
via the rp-process (Wallace & WoosleyQMl] Shatz et al. |2001|) . 
For these systems, Cumming & Bildsten (2001 ) have shown that 
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a small residual amount of 12 C (X I2 q * 0.1) in a heavy element 
bath might be enough to ignite a superburst. Indeed, the low ther- 
mal conductivity of the rp-process ashes gives a large tempera- 
ture gradient in the ocean which favours an unstable ignition of 
12 C (Cumming & Bildsten 2001) (By ocean, we mean the re- 
gion far below the zone where the accreted matter is decelerated 
from its free-fall velocity, and directly underneath the hydro- 
gen/helium burning layer (Brown & Bildsten [19981 0. Moreover 
Shatz et al. (2003) have shown that, at the high temperature 
reached in superbursts (T > 10 9 K), the photodesintegration of 
the heavy rp-process ashes releases a quantity of nuclear energy 
which can be larger than the energy release from the fusion of 
carbon. If this is true, superbursts are the only known cosmic 
phenomena where photodisintegration of heavy elements is the 
main energy source. This interesting nucleosynthesis aspect will 
be considered in a future work. 

This paper deals only with the case of pure He accretion. 
Pure He accretors are rare systems. There is only one super- 
burst which has been observed in a system where the accreted 
material likely has a very high He abundance. This system is 
4U 1820-30 (Strohmayer & Brown 120021 hereafter SB02). It is 
an ultra compact binary where the companion star is probably 
a low-mass helium dwarf (SB2002, Cumming 2003). This su- 
perburst detected by RXTE was * 3 hours long, had a peak lu- 
minosity of w 3.4 10 38 ergs s -1 , an observed energy release of 
1.5 10 42 ergs, and was preceded by a normal type I burst (20 s 
duration) (SB02). SB02 showed that this superburst is thermonu- 
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clear in origin, fueled by the burning of carbon produced by the 
stable burning of the accreted helium between bursts. As neutri- 
nos carry away most of the energy, the total energy of the super- 
burst must be much greater than the observed one (» 10 44 ergs). 
The ashes of He burning depend on the stability of the combus- 
tion. If stable, C is the main product (Brown & Bildsten 1998 ), 
while an unstable combustion during normal type I bursts pro- 
duces iron group elements (SB02). SB02 have shown that the 
iron made during bursts may mix with the carbon made during 
stable burning between bursts, so that the deep ocean is a mixture 
of the two. For this reason, we will study here the characteristics 
of a detonation wave in pure 12 C and in a mixture X i2 Q = 0.3 
and X 52 p e = 0.7, which are the limiting cases considered by 
SB02. Due to the lack of protons, no rp-process can develop, 
which implies the absence of photodisintegration of heavy rp nu- 
clides. In such conditions, a restricted nuclear reaction network 
is sufficient to describe the nucleosynthesis. This greatly eases 
the hydrodynamical simulations. We use a 13 species a chain 
network commonly used in other astrophysical hydrodynamics 
code (Fryxell et al. [19891 . 

Previous works concerning superbursts focused on the ther- 
modynamical state of the surface layers just before ignition, and 
on the cooling following the bursts (Weinberg et al. 2006), but 
detailed hydrodynamic calculations are still missing. The obser- 
vation of oscillations during the 2001 February 22 (UT) super- 
burst from 4U 1636 53 suggests some departure from spherical 
symmetry in the superburst phenomenon (Strohmayer & Brown 
2002b). Indeed, as it was suggested for normal type I X-ray 
bursts (Shara 1982 ), it is unlikely that ignition conditions will be 
achieved over the entire surface simultaneously. It appears more 
likely that burning is initiated locally, and then spreads laterally 
around the neutron star. Since the neutron star is rotating, the os- 
cillations might be understood with a model including rotation 
and a non-uniform surface brightness, as well as the spreading 
of the combustion around the surface (Spitkovsky et al. 120021 ). 



2. Numerical method 

Our simulations are performed with a modified version of the un- 
split, shock-capturing algorithm for multi-dimentional systems 
of hyperbolic conservations laws with source terms proposed by 
Papalexandis et al. (2002). It is a finite-volume method in the 
spirit of the original MUSCL scheme of van Leer (1979 ). The 
algorithm is of second-order in the smooth part of the flow. It 
avoids dimensional splitting. For the purpose of our astrophys- 
ical study, the original algorithm is extended to treat a stellar 
equation of state and a thermonuclear reaction network. We have 
implemented a Riemann solver based on the one of Colella & 
Glaz ( 1985 ), which is able to treat a general equation of state like 
an astrophysical one (Fryxell et al. 119891 ). The parallelization, 
deemed necessary in order to handle the computational require- 
ments for the problem in hand, is based on the " mpi " library, as 
described in Deledicque & Papalexandris (2006 ). 

2.1. Hydrodynamics 

The algorithm solves the adiabatic Euler's equations for com- 
pressible, non viscous gas dynamics with source terms in two 
dimentions. The equations can be written in conservative form 
as 

QT T 
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To close the system, we need an equation of state of the form 



A combustion front may propagate in different ways 
(Williams 1965). According to Weinberg et al. (2006) the su- 
perburst rise evolves through three nuclear burning stages: an 
hour-long convective stage, a runaway stage, and a hydrody- 
namic stage. A combustion wave forms and may propagate from 
the site of the runaway as a detonation. 

The only previous numerical studies of the propagation of 
a detonation front at the surface of a neutron star were made 
by Fryxell & Woosley ( [19821 ) and by Zingale et al. ( T200TT ). They 
considered only a detonation propagating in the He layer and not 
in the underlying C layer of relevance to superbursts. This paper 
presents the first hydrodynamical simulation of the C-detonation 
type. This simulation helps illustrate the performance of a new 
code we have developed. We limit ourselves here to the one- 
dimensional case. Two-dimensional results will be presented in 
a forthcoming paper. 

Our code is based on a new MUSCL-type parallelized al- 
gorithm introduced by Papalexandis et al. (2002) and extended 
to cope with astrophysical conditions. It is described in Sect. [2] 
Our time-dependent simulations of detonation in pure 12 C or in a 
mixture of 12 C and 52 Fe at constant pressure and density typical 
of superbursts are compared in Sect. [3] with steady state predic- 
tions. The importance of the resolution is discussed for different 
initial conditions. Section [4] contains our conclusions and per- 
spectives. 



p = p(p,T,Y), (3) 
e = e(p,T,Y). (4) 

In the equations above, p is the density, u = (u, v) is the velocity 
vector, p is the pressure, e, — e + y is the sum of the specific 
internal energy and the specific kinetic energy, F, is the molar 
fraction of species /Q, Y is the vector of F;, T is the temperature 
and s"" c is the total rate of thermonuclear energy released per 
gram of matter. For a general reaction as cj + CjJ ^ c^K + ciL, 
where cijxi are the stoechiometric coefficients of the nuclides /, 
J ,K and L, s m,c and fl™ c are defined by 

s m,c = -N A ^] M,c 2 /?™ e , (5) 
and 

rt = Yj a ^ AkYk+ Yj {ai(iJ)[i ' j]kY ^ Y 7 + 

k j,kj 

Oi(k, l)[k, q^Ff } + ai {k, k, k)p 2 N 2 A (crvh k Yl (6) 

1 Note that the molar fraction of the leptonic species Y L does not 
appear in equations 1314-1 In the thermodynamic conditions relevant to 
superbursts, the plasma can be considered as totally ionized. The con- 
dition of electroneutrality allows the calculation of Y L once the 7, 's are 
known. (Cox & Giuli [T9g8l ) 
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where [i, j] k = pNA{crv)i+j^>k> Na is the Avogadro number, n sp is 
the number of species in the system and M,c 2 is the rest mass 
energy of species i. The notation (crv),^^ represents the ther- 
monuclear reaction rate of the process i + j — » k + I per pair of 
particles The quantities a-, are statistical factors determined 
by the c iJX i (Arnett |l996l . 

2.2. The algorithm 

Equations ([TJ and (O take the integral form 



pu • dS = 0, 



dt f v pdV + Js 

— f pudV + f pun ■ dS+ f pdS = 0, 

* Jv Js Js 

— J pe,dV+ J pe,u dS + J pu-dS- J pe m,c dV = 0, 



4- f pYidV + f pYiu-dS 
dt Jv J s 



f p r T 

Jv 



dV = 0. 



(7) 
(8) 
(9) 
(10) 



These equations are written for an arbitrary control volume V 
whose boundary S has zero velocity. The hydrodynamical part 
of these equations is treated in the same way as in Papalexandris 
et al. (120021 1 . The nuclear part of the system requires a specific 
treatment, however. The original algorithm of Papalexandris et 
al. d2002) is indeed able to treat a single stiff source term. Instead 
our nuclear reaction network which comprises 27 reactions (see 
Sect. 12.3b introduces a set of stiff differential equations, and so 
very different time scales. This requires the adoption of a time- 
splitting version of the algorithm. We keep avoiding the dimen- 
sional splitting. The time splitting is of the Strang type (Strang 
1968). In this case the system of equations (fl]|2]i is solved by the 
split scheme 



tr 



(11) 



where U n+1 is the solution at time t+At. Here X/ is the numerical 
solution operator for the corresponding homogeneous conserva- 
tion law 



d_ 
dt 



U + V • F(U) = 0, 



(12) 



and Xj is the numerical solution operator for the system of ordi- 
nary differential equations 



dt 



U = G(U). 



(13) 



It is obtained from the semi-implicit extrapolation method of 
Bader & Deuflhard ( 119831 1 which is used to solve the nuclear 
part of the system of equations (fTlfTOli. 

The procedure of discretization and numerical evaluation of 
the hydrodynamical part of the integrals (TTTfTOb at each computa- 
tional cell is the same as in Papalexandris et al. (2002). However 
the gamma-law Riemann solver of the initial algorithm is re- 
placed by one based on the method of Colella & Glaz (1985) 
able to treat a general equation of state of the form given by 
equations ([3]|4|. 

The numerical scheme, which evaluates the solution at time 
(n + I) At from the solution at the previous time nAt for the hy- 
drodynamical part of the system (eq. IT2t can be written as 



(m u T 



(my)" - Af[(/n s • F m )^y 2 2 . - (ln s ■ F,„)^ .] 



s«+l/2 



At[Qn s • F m ) 



n+l/2 



y+i/2 



( /n 5 -F m )^ 2 ] 



(14) 



(mijUij)" + = (niijiiij)" - 

At[(ln s ■ V^jlj 
At[(ln s ■ F„)^g 2 • 



(15) 



AtWns ■ F v ) 



At[(ln s ■ F v ) 



n+l/2 

;+i/2,> ' 

n+l/2 
/J+l/2 ' 



(16) 



(mije tiJ ) 



(mtjetijT ~ 

Atms-F^ll-ilns-F^] 



(17) 



(m i jY iJ>k ) n+1 = (mijY idk r 
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where / and are the length of a cell interface and the unit 
vector normal to a cell interface, respectively. The flux vectors 
are given by 



F m = [pu,pv], 

F„ = [pu 2 + p,puv], 

F v = [puv,pv 2 + p], 

F e = [pe t u + pu,pe t v + pv], 

F Yi = [pY^pYiv]. 



(19) 
(20) 
(21) 
(22) 
(23) 



All the numerical simulations presented in this paper are per- 
formed with a Courant number (Leveque 1999) CFL = 0.3, and 
the temperature is not allowed to change by more than 10% dur- 
ing a time step (Fryxell et al. 1989). 

2.3. Equation of state and nuclear reaction network 

In the thermodynamic al conditions relevant to the surface of a 
neutron star, the plasma can be considered as fully ionized. Our 
equation of state accounts for partially degenerate and partially 
relativistic electrons and positrons. The ions are treated as a 
Maxwell-Boltzmann gas, and the radiation, considered to be at 
local thermodynamic equilibrium with the matter, follows the 
Planck law. We use a tabulated equation of state in the spirit of 
Timmes's one (Timmes & Arnett 1999). Coulomb interactions 
of the bare nuclei with the surrounding electron-positron gas 
are not taken into account here, and will be included in a future 
work. 

The selected nuclear reaction network is the one usually 
adopted in astrophysical hydrodynamics simulations in order to 
provide an energy source representative of explosive helium and 
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carbon burning in absence of hydrogen. It involves 13 nuclides 
( 4 He, 12 C, 16 0, 20 Ne, 24 Mg, 28 Si, 32 S, 36 Ar, 40 Ca, 44 Ti, 48 Cr, 52 Fe 
and 56 Ni) linked by 27 reactions comprising the 1 1 (a, y) re- 
actions from 12 C(ff, y) 16 to 52 Fe(a, y) 56 Ni, the correspond- 
ing 11 endothermic photodesintegrations, the three heavy-ion 
reactions 12 C( 12 C,a) 20 Ne, 12 C( 16 0,a) 24 Mg and 16 0( 16 0,a) 28 Si, 
and the triple-alpha reaction and its inverse. As in Fryxell et 
al. ([1989), the reaction rates are taken from Thielemann et al. 
dl9861 l, where each reaction rate is given in the temperature in- 
terval 10 8 < T < 10 10 K by 

N A (crv} = exp{c] + c 2 T 9 l + c 3 T 9 1/3 + c 4 T 9 l/3 + 

c 5 T 9 + c 6 T 5 9 /3 + c 7 ln(T 9 )), (24) 

where T 9 = 10 9 K and the values of the numerical coeficients q 
are given by the authors. The network equations are constructed 
as described by Eq. ( fT3l , and are solved with the use of the vari- 
able order Bader-Deuflhard semi-implicit time integrator (Bader 
& Deuflhard [19831 suggested by Timmes ( [19991 . The way this 
method is introduced in the algorithm is the same as in Press et 
al. (fl992l . 



i — 1 — 1 — 1 — i — • — • — • — i — <■ 




4 -2 D 2 + 



log(Z) 



2.4. Validation tests 

The algorithm has already been validated for terrestrial detona- 
tions with a gamma law equation of state, and a single chem- 
ical reaction kinetic represented by an simple Arrhenius law 
(Papalexandis et al. I2002.l t. Its validity in astrophysical situa- 
tions with a general equation of state and a nuclear reaction net- 
work has been checked with some validation tests considered 
by Fryxell et al. d 1989b . They involve non reactive and reactive 
shock tubes with the astrophysical equation of state and the nu- 
clear reaction network. Our results are similar to those of Fryxell 
et al. ( [19891 . 

The comparaison of our time-dependent results with a steady 
state calculation, as presented in Sect [3] also validate the accu- 
racy of our algorithm. 



3. Detonation profiles 

One-dimensional steady-state calculations (i.e. calculations 
where the detonation speed remains constant) provide the main 
parameters characterizing a detonation, such as the characteristic 
time and length-scales and the reaction-zone structure. These ba- 
sic detonation properties are necessary to set the initial parame- 
ters and boundary conditions in the time-dependent calculations. 
The treatment of the steady-state case is called the ZND model 
(Fickett & Davis 1979). According to this model, the detona- 
tion consists of an infinitely thin shock followed by a burning 
zone. All the reactions take place inside this zone and the re- 
leased energy sustains the shock. In the laboratoy frame, these 
steady-state equations are 



dp 

dt ~ (D- u) 2 



•7 



de _ f_dp +£nuc 
dt p 2 dt 



dZ p^ 
— = D - u — — , 
dt p 



(25) 

(26) 
(27) 



Fig.l. 

4 He, 12 C, 
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o 



Nuclear mass 
28 Si, 32 S, 52 Fe and 



fraction 
56 Ni for a 



pure 12 C at T = 10 8 K andp = 10 8 gem- 



profiles of 
detonation front 
Z is the distance 
to the shock in cm. The thin solid lines give the steady state 
solution of the ZND model. The heavy dots (the density of the 
dots is so high that they form thick curves) are obtained with 
the time-dependent calculations with two different resolutions 
(mesh size 10~ 4 cm and 10 cm). 



and 

dJl = Rnuc 
dt '' ' 



(28) 



where D, Z and po are the detonation speed, the distance behind 
the shock and the density of the unburnt matter. The quantity 

1 /2 

<t> = E mlc {dplde) p ,Y j is the thermicity, aj — (dp / 'dp) s r is the 
frozen sound speed and S is the entropy (Khoklov 1989). The 
system of equations ( 1251128b is closed by the equation of state 
(eqs. [3J3). For a given detonation speed, we compute a post- 
shock state using the Hugoniot relations (Khoklov l 19881 1 and we 
integrate from the shock to the end of the burning zone (when 
the R"" c vanish) using equations ( 12511281 ). 

Since superbursts in pure He accreting systems may be due 
to pure C detonation at high density and temperature we calcu- 
late the ZND profiles for a detonation in pure 12 C at a tempera- 
ture T - 10 8 K and a density p = 10 8 g cm 4 . The nuclear mass 
fraction profiles of some of the most abundant species are pre- 
sented in Fig.Q](thin solid lines). As can be seen a large variety 
of length scales are at work. The total reaction length given by 
the ZND model is of the order of 10 4 cm. Significant changes 
in the nuclear mass fractions occur already at 10~ 4 cm. This sug- 
gests that the full resolution of the detonation requires a time- 
dependent simulation over a domain as large as 10 4 cm with a 
resolution of 10~ 4 cm. Owe to computational time limitations we 
are unable to reach this resolution. However, following Gamezo 
et al. (119991 ), we have been able to perform two sets of calcula- 
tions, one with a resolution of 10 cm over a domain of 10 4 cm, 
and one with a resolution of 10 4 cm over a domain of 1 cm. 

The system of length / = 10 4 cm is considered with an inflow 
boundary condition at x = cm and an outflow boundary con- 
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dition at x =10 4 cm. As initial conditions, we select pure 12 C at 
T = 10 8 K, p — 10 8 g cirr 3 and a material velocity of cm 
s _1 . To trigger the detonation we set the initial conditions in an 
ignition zone between x — cm and x =10 3 cm to a temperature 
of 4.46 10 9 K, a density of 3.01 10 8 gcm -3 , a material velocity 
of 8.07 10 8 cms _1 and pure 56 Ni. We take 1000 numerical cells 
in the domain, leading to a resolution of 1cm. 

The simulation of a detonation of pure 12 C in the same ther- 
modynamic conditions as above has also been performed on a 
much smaller domain of lenght 1 = 1 cm with 1 4 numerical cells, 
so that a resolution of 10~ 4 cm is achieved. The initial disconti- 
nuity is positioned at 0. 1 cm. 

We have superimposed the profiles of the nuclear mass frac- 
tions of some of the most abundant species as a function of the 
distance to the shock obtained with the ZND algorithm and with 
the hydrodynamic algorithm. The ZND profiles are obtained for 
a velocity of propagation of the detonation wave D = 1.3 10 9 
cm s -1 , which is the velocity of the detonation front in the time- 
dependent hydrodynamic simulation. As can be seen in Fig. [T] 
the time-dependent profiles are very close to the steady state 
ones. So, as in Gamezo et al. ( 119991 ), the partial resolution ap- 
proach can be applied for the simulation of detonations in this 
system. This will be very useful in future 2D calculations where 
the computational time is crucial. 

The preceding comparison between the steady-state and the 
time-dependent results was made without taking nuclear reac- 
tion rate screening effects into account. However, at high densi- 
ties, screening may be important (Cox & Giuli 1968). and all the 
following results are obtained with the adoption of the screening 
corrections of Wallace et al. ( 1982). With this correction we have 
performed the simulation of a detonation of pure 12 C in the same 
thermodynamic conditions as above at four different resolutions: 
one with a domain of length 1 = 1cm, one with 1 = 100cm, one 
with 1 = 1000cm, and one with 1 = 10 4 cm, always using 1000 nu- 
merical cells. The four corresponding nuclear mass fraction pro- 
files are presented in Fig. [2] One sees the impact of the screening 
effects by comparing Figs. Q] and [2] Screening decreases the re- 
action length, increases the detonation velocity, and modifies the 
final composition. However these effects are small (less than 1 % 
of the detonation velocity). 

The profiles of temperature, velocity, density and pressure at 
time t = 5 10" 6 s for a resolution of 10 cm are presented in Fig. 
[3] The same profiles, and the nuclear energy generation profiles 
at time t = 6 10~ 10 s for a resolution of 0.1 cm are presented in 
Figs. |4] and [5] Ninety percent of the nuclear energy is already 
liberated over a distance of 10~ 2 cm. 

If we want to study the propagation of the detonation around 
a neutron star surface, we will be mainly interested in the struc- 
ture of the end of the reaction zone. Indeed the detonation travels 
approximately 6 10 6 cm (R^s ~ 10 km). What happens on 10~ 2 
cm is thus quite irrelevant, even if 90% of the nuclear energy is 
alredy liberated there. With a resolution of only 10 cm, we can 
already infer the global properties of the detonation in spite of 
the fact that we miss the destruction of C, the production and 
destruction of O, the production of Si and of S. All these reac- 
tions occur within one computational cell just behind the shock, 
but this has no global impact on the mean thermodynamic value 
at the end of the reaction length. This can be seen by comparing 
the thermodynamic profiles obtained with both resolutions (Figs. 
[3]and[4](. However, fully resolving the detonation would be very 
important for the study of the extinction of the detonation (Maier 
& Niemeyer 2006). This question is not tackled here. The pro- 
files of Fig. Ogive approximate values only, as the detonation is 
not perfectly resolved. 




Fig. 2. Same as Fig.Q] but with screening effects taken into ac- 
count. Results only obtained with the time-dependent calcula- 
tions for four resolutions (from left to right: mesh size 10 3 cm, 
0.1 cm, 1 cm and 10 cm). 



A simulation with a composition X n Q = 0.3 and X 5: p e = 0.7 
(SB02) has also been conducted, the ignition conditions being 
the same as for pure 12 C. The temperature, velocity, density and 
pressure profiles at time t = 7 10~ 6 s are presented in Fig. [6] In 
Fig. Q the nuclear mass fraction of some species is shown as a 
function of the distance to the shock. We have only performed 
a simulation with a 10 cm resolution for the reasons explained 
above. 

From these two sets of results (detonation in pure 12 C and in 
a mixture of 12 C and 52 Fe), we see that the composition of the 
material before the passage of the detonation wave affects the 
velocity of propagation of the detonation wave and the composi- 
tion in the burned material. For pure 12 C the velocity of propaga- 
tion of the detonation is D » 1.3 10 9 cm s _1 , and for the mixture 
of 12 C and 52 Fe, D * 1.22 10 9 cm s" 1 . The pure 12 C detonation 
mainly produces 4 He, and the mixed detonation leads essentially 
to 56 Ni. 



4. Conclusions 

Our hydrodynamical algorithm for modeling astrophysical det- 
onations has been shown to be robust to test cases. In partic- 
ular it reproduces quite well the steady state solution obtained 
with a totally different code. This give us confidence for future 
multi-dimensional simulations. Some improvements still need to 
be made, like the inclusion of gravity, of non-ideal terms in the 
equation of state, and more up-to-date data for the nuclear reac- 
tion network. 

We have underlined the large difference between the total 
reaction length and the length on which some species (e.g. car- 
bon) burn in conditions relevant to superbursts. This difference 
leads to enormous numerical difficulties because all the length 
scales cannot be resolved at the same time during a single sim- 
ulation (except possibly with some sub-grid models). We have 
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Fig. 3. Temperature (in K), velocity (in cms -1 ), density (in 
g cm 4 ) and pressure (in erg crrT 3 ) profiles of a detonation front 
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Fig. 5. Nuclear energy generation (ergg -1 s~') profile of a deto- 
nation front in pure 12 C at T — 10 8 K and p - 10 8 g cm 4 at time 
= 6 10~ 10 s with a resolution of 10~ 3 cm. X is in cm. 
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shown that the carbon detonation in superburst conditions might 
be studied by a partial resolution approach. This important, es- 
pecially in multidimensional simulations where the computation 
time is a crucial limitation. 

Our simulations give the global thermodynamic state of the 
material after the passage of the detonation in the carbon layer 
at the surface of a neutron star. These conditions are not very 
sensitive to the exact initial composition of the matter, in contrast 



to the final composition. In both cases, however, all the carbon 
is burned. It has to be replenished to allow the occurence of a 
subsequent superburst. 

It is unlikely that the entire carbon layer ignites at the same 
time. More probably, ignition spots develop. Since the detona- 
tion velocity depends on the composition of the carbon layer 
before ignition, the time for the detonation to propagate laterally 
around the neutron star also depends on the composition. 
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In a subsequent paper we will investigate a larger parameter 
space of thermodynamic condition. In particular, we will study 
the impact of the initial temperature or density. It would also 
be more realistic to introduce an initial temperature and density 
profile. Two-dimentional superburst simulations are under study 
and will complement the work of Zingale et al. ( 2001) for normal 
type I X-ray bursts. We will also investigate the vertical prop- 
agation of the detonation and its interaction with an overlying 
helium layer. One of the aims of this simulation is to examine if 
the penetration of a detonation wave into the helium layer could 
given rise to the precursor observed prior to the superburst from 
4U 1820-30. 
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at time = 7 10 s. Z is the distance to the shock in 



